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This paper summarises a numerical investigation of the statistical properties of orbits evolved in 
'frozen,' time-independent iV-body realisations of smooth, time-independent density distributions 
corresponding to integrable potentials, allowing for 10 2 ' 5 < N < 10 5 ' 5 . Two principal conclusions 
were reached: (1) In the limit of a nearly 'unsoftened' two-body kernel, i.e., V(r) oc (r 2 +e 2 ) -1 ^ 2 for 
small e, the value of the largest Lyapunov exponent x does not appear to decrease systematically 
with increasing N, so that, viewed in terms of the sensitivity of individual orbits to small changes 
in initial conditions, there is no sense in which chaos 'turns off' for large N. (2) Nevertheless, 
there is a clear, quantifiable sense in which, on the average, as N increases chaotic orbits in the 
frozen- N systems come to more closely resemble integrable characteristics in the smooth potential. 
When viewed in configuration or velocity space, or as probed by collisionless invariants like angular 
momentum, frozen- N orbits typically diverge from smooth potential characteristics as a power law 
in time on a time scale oc N p to, with to a, characteristic dynamical, or crossing, time. For the case of 
angular momentum, the divergence is well approximated by a t 1 / 2 dependence, so that, when viewed 
in terms of collisionless invariants, discreteness effects acts as a diffusion process which, presumably, 
can be modeled by nearly white Gaussian noise in the context of a Langevin or Fokker-Planck 
description. For position and velocity, the divergence is somewhat more rapid and characterised by 
a t q power law growth with q « 1, a result that likely reflects the effects of linear phase mixing. 



PACS number(s): 05.60.+W, 51.10.+y, 05.40.+j 
I. INTRODUCTION AND MOTIVATION 

Many astronomical objects, including, e.g., globular 
clusters, are typically modeled by bulk gravitational po- 
tentials which manifest a high degree of symmetry and 
which, being integrable, lead to completely regular char- 
acteristics with no possibility of chaotic behaviour. One 
knows, however, that such bulk potentials constitute ide- 
alisations, the true system corresponding (at least ap- 
proximately) to a realisation of the gravitational iV-body 
problem. The important point, then, is that motion in 
the iV-body problem, even for an iV-body system which 
samples a smooth, time-independent phase space distri- 
bution corresponding to an integrable potential, is typi- 
cally chaotic in the sense that orbits exhibit exponential 
sensitivity towards small changes in initial conditions jj] . 
This perhaps is not surprising. The true potential as- 
sociated with a collection of point masses no longer pos- 
sesses the symmetries of the original integrable potential, 
so that there is no reason why the orbits should not be 
chaotic. 

However, what is, perhaps, surprising is the expecta- 
tion, derived both from theoretical arguments ^|J|] and 
from numerical simulations B , that the TV-body problem 
remains chaotic even for very large N. Suppose, e.g., that 
a system of total mass M = 1 is represented by a collec- 
tion of N objects of mass m — l/N, so distributed as 
to sample the density distribution corresponding to an 



integrable potential. The claim then is that, when ex- 
pressed in units of a natural dynamical, or crossing, time 
t d ~ 1/ \[Gp, with p a typical density, the characteristic 
time scale t on which an initial perturbation in any given 
orbit tends to grow will not diverge for N — > oo. In this 
sense, the degree of chaos manifested by individual orbits 
is not expected to 'turn off' for very large N. There is an 
apparent consensus, motivated both from theory and nu- 
merical experiments, that r should not increase without 
bound for N — *■ oo, although there is some disagreement 
in the literature as to whether t(N) should converge to- 
wards an iV-independent value |j or whether r should 
instead slowly decrease with increasing N |J . 

If, however, this be true, one is confronted with subtle 
questions of principle regarding the nature of the contin- 
uum limit. It is generally assumed [9 that, for sufficiently 
large N, a self-gravitating system of discrete point masses 
can be characterised adequately by a smooth phase space 
density that solves the collisionless Boltzmann equation 
(CBE), i.e., the gravitational analogue of the Vlasov 
equation from plasma physics. The obvious point, then, 
is that time-independent solutions to this equation which 
manifest a high degree of symmetry correspond typically 
to bulk potentials which are integrable or, even if they 
be nonintegrable, admit large measures of regular orbits. 
But how is one to reconcile integrable or near-integrable 
behaviour in such bulk potentials with the presumed fact 
that, even for very large N, individual orbits in the true 
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iV-body problem typically manifest chaotic behaviour on 
a time scale ~ £d? 

Strictly speaking, there is no logical contradiction: It 
is completely possible for collective properties of an N- 
body system to be described correctly by the CBE, even 
if the characteristics associated with the self-consistent 
potential do not coincide, even approximately, with real 
iV-body trajectories. However, it would seem important 
to pin down carefully what is actually going on: 

• Is it really true that individual trajectories in the N- 
body problem are chaotic for very large N, even if the 
bulk potential associated with the system is integrable? 
The indications are that the answer to this is: yes. How- 
ever, most of the work done to date on chaos in the N- 
body problem has focused on systems with comparatively 
small N and/or a hierarchy of masses, or, for larger sys- 
tems, on comparatively short time behaviour. Little if 
any work has been done to provide estimates of honest 
Lyapunov exponents over intervals 3> tjy for large N sys- 
tems comprised of bodies of comparable mass. 

• Does this chaos reflect simply the fact that the sys- 
tem is grainy, or does it reflect the details of the full 
iV-body dynamics? If, e.g., one were to replace a smooth 
time-independent, integrable potential by the irregular 
potential associated with an iV-particle sampling of the 
smooth mass density that is frozen in time, to what ex- 
tent will motion in that grainy time-independent poten- 
tial be chaotic? 

• Even presuming that the TV-body problem is chaotic 
on a time scale *~ tjj, is there some well defined sense 
in which the 'average' properties of individual iV-body 
trajectories track characteristics given by the CBE? 

These conceptual issues are also related directly to the 
problem of 'softening.' It is generally recognised that, 
for small TV, close encounters between individual masses 
are more important dynamically than for larger TV ||. 
For this reason, TV-body simulators interested in explor- 
ing the physics of the TV-body problem for larger TV often 
suppress the effects of close encounters artificially by re- 
placing the true 1/r potential by a softened potential 
V(r)(x (r 2 + e 2 ) -1 / 2 for some 'softening parameter' e. 
This certainly suppresses encounters with impact param- 
eters < e which, presumably, is a good thing. However, 
there are strong indications that orbits integrated with 
such a softened potential tend to be 'less chaotic' in their 
behaviour, so that the introduction of softening also has 
the potentially undesirable effect of removing TV-body 
chaos which really ought to be present, even for very 
large TV. In any event, earlier investigations of chaos in 
the TV-body problem based on simulations which incor- 
porate a large amount of softening must be viewed with 
suspicion, since such simulations could suppress precisely 
the effects which one might wish to explore! 

This paper summarises a detailed exploration of 
chaos in time-independent potentials generated by sam- 
pling the smooth density p(r) associated with a time- 
independent solution to the CBE to create a frozen TV- 
body realisation of that equilibrium. Most of the work 



focuses on the particularly simple case of an integrable 
Plummer potential which derives from a spherically 
symmetric mass distribution. However, it was also con- 
firmed that, modulo one point discussed in the conclud- 
ing section, the same qualitative results obtain for the po- 
tential associated with a constant density spherical con- 
figuration. 

Section II begins by describing the numerical exper- 
iments that were performed. Section III summarises a 
computation of honest Lyapunov exponents in frozen TV- 
body realisations of the Plummer potential, exploring 
how the largest exponent x associated with represen- 
tative initial conditions varies as a function of e and 
TV. The principal conclusion here is that, at least for 
small values of e, orbits in such potentials are invariably 
chaotic; and that, even for particle number as large as 
N = 10 5 ' 5 , there is no sense in which increasing TV 'turns 
the chaos off.' Section IV demonstrates that, even though 
the Lyapunov exponents do not decrease with increasing 
TV, there is a well defined sense in which, as TV increases, 
orbits in frozen-TV potentials remain 'close' to smooth 
potential characteristics with the same initial condition 
for progressively longer times. Section V concludes by 
summarising the principal conclusions, providing a sim- 
ple physical interpretation, and then commenting on po- 
tential implications. 

The principal conclusion of this paper is that, for inte- 
grable smooth potentials which admit no chaos, the con- 
tinuum limit makes sense even at the level of pointwise 
properties of individual trajectories. The possibility of 
chaotic characteristics leads necessarily to very different 
behaviour and, for this reason, the case of nonintegrable 
potentials that admit both regular and chaotic charac- 
teristics will be considered in a separate paper. 



II. DESCRIPTION OF THE NUMERICAL 
EXPERIMENTS 

The numerical computations reported here were per- 
formed for a so-called Plummer potential, 
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This potential is generated via Poisson's equation from a 
density profile 
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and corresponds to an equilibrium solution to the CBE 
satisfying 

/ A{-E) 7 / 2 if *(r = 0) < E = \v 2 + $ < 0; 
H j_ \ if £=±« 2 + $>0. 

(2.3) 
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Units were so chosen that G = M = b = 1. 

The principal aim was to compare orbits generated 
in the smooth potential with orbits evolved in time- 
independent iV-body realisations of the potential. For 
a variety of fixed values of N and e, 20 different time- 
independent iV-body potentials were constructed. Each 
of these was associated with a random sampling of the 
smooth density distribution generated using a von Neu- 
mann rejection algorithm (cf. Q). This entailed con- 
structing singular density distributions 
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which, allowing for a softening parameter e, yielded po- 
tentials of the form 
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The objective then was to select individual initial con- 
ditions (i"o,vo) and to evolve these initial conditions in 
both the smooth potential and the 20 'frozen' TV-body 
potentials, while simultaneously tracking the evolution 
of a small initial perturbation, periodically renormalised 
at fixed intervals St, so as to extract an estimate of the 
largest (short time) Lyapunov exponent || . 

The integrations were performed for a time correspond- 
ing physically to ~ lOOtn using a Runge-Kutta integra- 
tor that typically conserved energy to at least one part 
in 10 4 . The value lOOtrj was selected (i) because it corre- 
sponded to an interval sufficiently long that one began to 
see convergence towards a well-defined Lyapunov expo- 
nent x and, perhaps more importantly, (ii) because, for 
physical systems like real galaxies, lOOto corresponds to 
an interval comparable to the age of the Universe. 

The total particle number N in the 'frozen' TV-body 
potentials was allowed to vary between N = 10 2 ' 5 and 
N = lO 5 5 . Physical interest focuses primarily on the 
limit e — > 0, this corresponding to an 'honest' iV-body 
calculation. However, the effects of a nonzero e were 
also considered in some detail, with the aim of ascertain- 
ing possible undesirable consequences for conventional 
iV-body simulations, which typically involve a substan- 
tial softening. The experiments with variable e indicated 
that, for e < 10~ 4 or so, the precise value of e was largely 
immaterial, at least statistically. 



III. SHORT TIME LYAPUNOV EXPONENTS 

The principal diagnostic here was the mean (short 
time) Lyapunov exponent (x), generated, for a given 
choice of initial condition and for specified values of e and 
N, as the average value of x at t = lOOtn for 20 different 
frozen-iV potentials. The fundamental question was how, 
for fixed initial condition, this (x) depends on e and N. 



FIGURE 1 exhibits (x) as a function of log 10 e for mul- 
tiple integrations of one representative initial condition, 
which corresponded initially to a roughly isotropic distri- 
bution of velocities, allowing for several different values 
of N. FIGURE 2 gives (x) as a function of log 10 N for the 
same initial condition, now allowing for several different 
values of e. 

It is evident from FIGURE 1 that, at least for compar- 
atively large values of softening parameter, decreasing e 
tends to make the orbit more chaotic. This is hardly 
surprising: Since the smooth potential is integrable, one 
anticipates that the chaos is associated completely with 
close encounters between the test mass and individual 
frozen masses. The introduction of a nonzero smoothing 
corresponds de facto to the introduction of a minimum 
impact parameter (since the potential is bounded in mag- 
nitude by V max = —GM/Ne) but the existence of such 
a minimum impact parameter limits the maximum effect 
that can arise from a close encounter. 

However, for sufficiently small values of e, the precise 
value of e appears to be largely immaterial. This again 
is hardly surprising: As long as e is small compared with 
the value of the closest separation between the test parti- 
cle and any of the frozen particles during the course of the 
integration, the test particle feels an essentially unsoft- 
ened potential and should behave (at least statistically) 
as if e = 0. The point then is that, for N < 10 6 and an 
integration time as short as lOOt^, the minimum sepa- 
ration associated with the closest encounter between the 
test mass and any of the frozen masses should be greater 
than or comparable to e~ 10 . Indeed, a simple geo- 
metric argument indicates |J that the time scale t t for a 
close encounter with minimum separation as small as e 
scales as 
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(3.1) 



where R sys is the size of the system in question. 

One obvious implication of these results is that the in- 
troduction of a large amount of softening into a numerical 
simulation can have the unnatural result of significantly 
decreasing the amount of chaos manifested by individual 
orbits in a real astronomical system. 

For comparatively large values of e, (x) decreases 
rapidly with increasing TV but, for sufficiently small val- 
ues of N, it appears that (x) is nearly independent of 
N (although there are hints that (x) may continue to 
increase very slowly). The fact that, for large e, (x) 
should decrease with increasing N can again be explained 
by comparing the magnitude of e with the typical dis- 
tance between masses in the system, which is of order 
n^ 1 / 3 ~ Rsyg/N 1 / 3 , with n a characteristic number den- 
sity. If e is larger than, or comparable to, n" 1 / 3 , even 
weak close encounters are essentially 'turned off,' so that 
the source of chaos has been largely reduced, if not com- 
pletely removed. The fact that (x) should be essentially 
independent of N in the limit e — > has been argued 
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by various authors in a number of different ways . A 
simple physical explanation is provided in the concluding 
section. 

If a single orbit be integrated for progressively longer 
times, how quickly will the short time Lyapunov expo- 
nent x(t) converge towards the true time-independent x? 
Studies of orbits in smooth nonintegrable potentials re- 
veal that, when the phase space is highly complex and, 
because of the Arnold web, orbits can be 'stuck' tem- 
porarily in regions where the short time Lyapunov ex- 
ponents are especially small or especially large, the time 
required for a reasonable level of convergence can be ex- 
tremely long, ~ I0 5 t]j — 10 6 t]j or even more ||. If, how- 
ever, the phase space is simpler in the sense that the 
Arnold web forms less of an impediment and such trap- 
ping is comparatively rare, the time required is typically 
much shorter. One way in which to quantify the overall 
rate of convergence is by performing a simple time se- 
ries analysis: An orbit segment of length T can of course 
be divided into k segments of length At = T/k and a 
short time Lyapunov exponent x(Ai) computed for each 
segment. The dispersion a x (At) then provides a useful 
probe of the degree to which, on time scales ~ At, the 
degree of chaos exhibited by different orbit segments is 
more or less the same. Determining o~ x as a function of 
At provides a quantitative characterisation of the rate of 
convergence towards a unique Xoo- A simple argument 
based on the Central Limits Theorem suggests JkJ that, 
if the accessible phase space regions are simple and trap- 
ping is rare, so that the amounts of chaos manifested at 
times t and t + At are essentially uncorrelated, 

a x cx (At)~P, (3.2) 

with p — 1/2. If, alternatively, the phase space is com- 
plex and trapping is important, one would expect that 
a x decreases much more slowly with increasing At. 

Such a time series analysis was performed for the data 
sets used to generate the mean exponents (%). For 
each set of 20 integrations, each orbit segment of length 
T = 10(Md was separated into k segments of length 
At = T/k. A short time Lyapunov exponent x(At) 
was then computed for each of the resulting 20k seg- 
ments, and these were used to compute the dispersion 
<J x {At). Allowing for k = 2«, for q = 0, 1, 2, 3, 4, 5, and 6 
was equivalent to varying At between At = (100/64)to 
and At — lOOtz?. This time series analysis led to the 
conclusion that the dispersion a x is typically well fit 
by a power law dependence of the form given by eq. 
(3.2), although the exponent p tends to be somewhat 
smaller than p = 1/2, the best fit value typically satisfy- 
ing p ~ 0.4. Several examples are exhibited in FIGURE 
3. The fact that p is comparatively close to 1/2, rather 
than the much smaller values that are often observed in 
very 'sticky' nonintegrable potentials corroborates 
the intuition that, because the chaos in this problem is 
associated exclusively with close encounters, trapping is 
rare and the degree of chaos exhibited at different times 
tends to be statistically uncorrelated. 



IV. COMPARISON OF SMOOTH AND N-BODY 
ORBITS 

It is clear that, for sufficiently short times, a frozen- 
N orbit will coincide almost exactly with the smooth 
potential characteristic associated with the same initial 
condition. And similarly, it is clear that, at sufficiently 
late times, the irregularities in the frozen- A potential will 
cause the frozen-A orbit to deviate significantly from the 
smooth characteristic. Probing the validity of the con- 
tinuum limit at the level of individual orbits thus de- 
volves into determining the rate at which the frozen-A 
orbits and smooth characteristics diverge. In this con- 
nection, two obvious questions arise. Do frozen-A orbits 
diverge from the smooth characteristics exponentially or 
as a power law in time? And how does the overall rate 
of divergence depend on A? 

Such probes of the validity of the continuum limit dif- 
fer from the ordinary point of view, where convergence 
is typically defined in terms of quantities like bulk mo- 
ments of the system, ignoring completely the behaviour 
of individual trajectories. A possible intermediate char- 
acterisation is to focus not on the pointwise behaviour 
of the chaotic orbits but, instead, on quantities which 
might be less sensitive to the A-body chaos. In partic- 
ular, one can also ask: How do frozen-A orbits deviate 
from smooth characteristics in terms of quantities which, 
in the smooth potential, correspond to time-independent 
constants of the motion, like angular momentum in a 
spherically symmetric system? 

These questions were addressed here both visually and 
quantitatively through a computation of the statistical 
properties of frozen-A orbits. Given n = 20 different 
trajectories {{ri(t), Vj(f))}, i = l,...,n, and the smooth 
characteristic (r s (t),v s (t)) associated with the same ini- 
tial condition, there are two types of moments which one 
might choose to consider. Quantities like 

TCI 

(r> = -][> (4-1) 

»=l 

and 

Dt* = {\ ri -{T)\ 3 ) = ±^2\ ri -(T)\* (4.2) 

i 

and the corresponding quantities generated from v and 
J = rxv focus on the frozen-A orbits in and of them- 
selves. Alternatively, such moments as 

Ar 2 ^(|r-r s | 2 ) = i£|r,-r s | 2 (4.3) 

i 

and 

fr 2 = (|(r)-r s | 2 ) (4.4) 

compare the frozen-A orbits with the smooth potential 
characteristic and, as such, their behaviour as a function 
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of N is particularly relevant in understanding the con- 
tinuum limit. Overall, the quantities Dr, Sr, and Ar 
were found to exhibit comparatively similar evolutions, 
so that attention below focuses on the moments (r) and 
Ar, which seem especially natural physically. 

The most striking conclusion is that individual frozen- 
TV orbits typically diverge from the smooth characteris- 
tic as a power law in time, not exponentially. This is 
true both for comparatively large values of e, where the 
frozen- TV orbits are nearly regular, and for smaller values 
of e, where the orbits are much more chaotic. This result 
is perhaps surprising. One might naively have supposed 
that, since the frozen-iV orbits are strongly chaotic, at 
least for small e, they would tend to diverge exponentially 
from the smooth characteristics on a time scale r ~ X 1 - 
However, such an exponential divergence is most defi- 
nitely not observed. 

For how long does this power law divergence persist? 
Does it cease when the distance between the frozen-TV 
orbit and the smooth characteristic is still small, or does 
the divergence continue until the frozen- TV orbit and the 
smooth characteristic tend to be widely separated in con- 
figuration space? If, e.g., this divergence terminated at 
comparatively small separations, much smaller than the 
size of the system, one could argue that, even though the 
frozen-TV orbits are chaotic, they still remain 'close' to 
the smooth characteristics. The answer here is that this 
divergence continues until the typical separation has be- 
come comparable to the size of the configuration space 
region to which the orbits are confined and the frozen- iV 
orbit has become completely 'decorrelated' in appearance 
from the smooth potential characteristic. 

The same conclusion also obtains if one focuses on an 
ensemble of 20 frozen- TV orbits and probes their statisti- 
cal properties. The six panels of FIGURE 4, generated 
for the initial condition used in FIGURES 1-3, com- 
pare (x) for frozen- TV ensembles with the smooth x s for 
orbits evolved with e = 10 , allowing for six values of 
TV extending from TV = 316 to TV = 100000. In each 
case, one finds that, for sufficiently large t, (x) — > 0, as 
would be expected if the frozen- TV orbits have become 
completely different from one another and move through 
configuration space with random orientations. FIGURE 
5 compares the radial coordinates (r) and r s for the ex- 
treme case of an initial condition corresponding in the 
smooth potential to a purely radial orbit. 

The time scale to on which the frozen- TV orbits diverge 
from the smooth characteristic, and hence the time scale 
on which Ar grows, increases with increasing TV. Even 
though the frozen-TV orbits remain 'equally chaotic' in 
the sense that their Lyapunov exponents \ remain nearly 
constant, they remain close to the smooth characteristic 
for progressively longer times. 

The left hand panels of FIGURE 6 exhibit Ar/2 1 / 2 R S , 
with R 2 the mean value of r 2 associated with the smooth 
characteristic, as computed for the same initial condition 
evolved with e = 10~ 4 for TV = 1000 and TV = 100000. 
The right hand side exhibits the same data, recorded 



at intervals of 0.025irj, once they have been subjected 
to a boxcar averaging over an interval St = Ud- The 
large envelops associated with the curves in the left hand 
panels reflect, e.g., the fact that, at late times, individual 
orbits in the n — 20 orbit ensembles are oscillating about 
a value of unity. 

That Ar/2 1 ' 2 R S converges towards unity is a reflection 
of the fact that the orbits have indeed become completely 
different one from another: Given that the frozen- TV or- 
bits conserve energy, one might expect that their average 
distance from the origin should, on the average, be the 
same as for the smooth characteristic, so that 



(r 2 (t))^R 2 for t^oo. 
Assuming, however, that this be true and that 



<r(t)T.(t)) - 



for 



t — > oo, 



(4.5) 



(4.6) 



one infers that Ar 2 — > 2R 2 . Analogous behaviour is ob- 
served for the quantity Av/2 1 ' 2 V S , with V 2 defined cor- 
respondingly for the smooth characteristic. 

As is manifested by FIGURE 6, the growth of Ar and 
Av is roughly linear in time. Indeed, when comparing an 
ensemble of frozen-TV orbits with a smooth orbit char- 
acteristic generated from the same initial condition, one 
finds that, for small e, Ar and At; are both reasonably 
well fit by a linear growth law of the form 



Ar 



Av _ t 

V s to' 



(4.7) 



The growth time to, which is the same for both Ar and 
Av, satisfies 



t G ^A G Nn D , 



(4- 



with Ac of order unity and p ss 1/2. FIGURE 7 exhibits 
log 10 (tG/*u) as a function of log 10 TV for two different 
initial conditions evolved with e = 10 -5 . 

The fact that to scales as TV 1 / 2 would suggest that the 
divergence of the frozen- TV orbits from smooth character- 
istics reflects a diffusion process, associated with a collec- 
tion of random close encounters. However, this might in 
turn suggest that Ar and At; should grow as t x l 2 , rather 
than the approximately linear growth that was observed 
in the numerical simulations. Interesting, though, such a 
t 1 / 2 behaviour does obtain for quantities like angular mo- 
mentum, which are conserved absolutely in the smooth 
potential. Indeed, one finds that, for small e, A J satisfies 



AJ 2 

IT 



t 



(4.9) 



where J s is the typical magnitude of the angular momen- 
tum associated with a characteristic with the specified 
energy. The growth time tj again scales as TV 1 / 2 , but 
tends to be somewhat larger than to, so that 



tjRiAjN*t D , 



(4.10) 
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with Aj~3A G and p sal/2. FIGURE 8 exhibits 
log 10 (tj/tu) as a function of log 10 TV for the same in- 
tegrations used to generate FIGURE 7. 

There is also a clear visual sense in which, as TV in- 
creases, the frozen-TV orbits become progressively more 
regular in appearance. This is, e.g., evident in FIGURE 
9, which exhibits the x-y projections of representative 
frozen-TV orbits with TV varying between TV = 316 and 
TV = 316228, all generated from the same initial condi- 
tion and integrated for a time t = 25tp with e = 1CU 5 . 
The final panel exhibits the smooth characteristic asso- 
ciated with the same initial condition. The most obvious 
point is that, as TV increases, the configuration space re- 
gion to which the orbit is restricted more closely coincides 
with the region occupied by the characteristic. For ex- 
ample, only for the three largest values of TV is the orbit 
'centrophobic' in the same sense as the characteristic. 

Also evident is the fact that the orbit 'looks smoother' 
for larger values of TV. This visual impression reflects 
the fact that, as TV increases, the power associated with 
the Fourier spectrum of an orbit tends to become more 
concentrated near a few special frequencies. (Since the 
smooth orbit associated with the same initial condition 
is regular, all its power is concentrated at a countable set 
of discrete frequencies.) This trend is illustrated in the 
eight panels of FIGURE 10, each of which exhibits |a;(o;)| 
for a single frozen- TV orbit generated from the same initial 
condition. In each case, the data are so normalised that 
the peak frequency has |x(o;)| = 1. The spectra were 
generated from a time series of 4001 points, recorded at 
intervals of 0.025<£>. 

The degree to which the orbits become more nearly 
regular with increasing TV can be quantified by determin- 
ing |]ll"f the 'complexity' of the orbits, i.e., the number of 
frequencies in the discrete Fourier spectrum which con- 
tain an appreciable amount of power. Two such measures 
of complexity are illustrated in FIGURE 11, which was 
computed for ensembles of frozen- TV orbits with varying 
TV, all evolved with e = 10~ 5 and generated from the 
same initial condition. The solid curve exhibits /o.i, de- 
fined as the sum of the numbers of frequencies fo.i.x, 
fo.i,y, and fo.iz which have more than 10% as much 
power as the peak frequencies for u> x , lu v , and u> z , i.e., 

fo.i = fo.i,x + fo.i,y + /o.i.z- (4-11) 

The dashed curve curve exhibits fco.95, defined corre- 
spondingly as the sum of the numbers of frequencies re- 
quired to capture 95% of the power in the x-, y-, and z- 
directions. In each case, the curve represents an average 
over different orbits in the ensemble, and the error bars 
represent the associated dispersions. The obvious point 
is that both these quantities decrease with increasing TV. 

The fact that, as TV increases, power becomes more 
concentrated near a few special frequencies has important 
implications for various physical processes which rely on 
resonances. For example, a variety of recent arguments 
in both galactic and solar system dynamics invoke a pro- 
cess of so-called of 'resonant relaxation,' p2] which relies 



on the assumption that, in the presence of a large cen- 
tral object (a supermassive black hole in the center of 
a galaxy or the Sun at the center of the solar system), 
TV-body orbits behave very nearly as if they were Keple- 
rian trajectories in the fixed 1/r potential associated with 
the central object. If the chaos exhibited by individual 
orbits jl3| implied that these orbits were highly irregu- 
lar, so that their power was not concentrated near the 
special Keplerian frequencies, resonant relaxation might 
seem quite implausible. Given, however, that the or- 
bits become progressively more regular for increasing TV, 
resonant relaxation would seem eminently reasonable, at 
least for systems in which TV is sufficiently large. 

V. CONCLUSIONS AND DISCUSSION 

Even though trajectories remain chaotic in the sense 
that the largest Lyapunov exponent does not decrease 
towards zero, there is a clear sense in which, for increas- 
ing TV, orbits in frozen- TV potentials exhibit a pointwisc 
convergence towards characteristics in the smooth poten- 
tial which the frozen- TV potentials were chosen to sample. 
Viewed in configuration or velocity space, frozen- TV orbits 
tend to diverge linearly from the smooth characteristic 
with the same initial condition on a time scale t G that is 
proportional to TV 1 / 2 . In this sense, the continuum limit 
appears justified even at the level of individual trajec- 
tories. The fact that frozen- TV orbits remain chaotic for 
very large TV is completely consistent with the existence 
of a well-defined continuum limit. 

It is easy to understand qualitatively why the frozen- 
TV orbits should remain chaotic even for very large TV. 
Given that the chaos disappears completely in the con- 
tinuum limit, where the orbits reduce to integrable char- 
acteristics, it would seem clear that the chaos must be 
associated with a sequence of 'random' interactions be- 
tween a 'test' particle and a collection of 'field' particles. 
However, this would suggest that the time scale associ- 
ated with the growth of a small initial perturbation can 
be estimated by considering the tidal effects associated 
with a pair of particles separated by a distance compa- 
rable to the typical interparticle separation. This tidal 
acceleration will of course scale as 

Sr = (<Sr-V)a ~ -^5r, (5.1) 

with r the separation and m the particle mass. Given, 
however, that r ~ n -1 / 3 ~ TV -1 ' 3 .^^, with n a charac- 
teristic number density and R sys the size of the system, 
it follows that the time scale t* associated with the in- 
teraction should satisfy 

U ~ 1/y/Gp. (5.2) 

In other words, the time scale associated with any orbital 
instability induced by the graininess of the system should 
be comparable to the dynamical time to, independent 
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of particle number TV. As TV increases, the size of the 
individual particle mass, m, and the cube of the typical 
separation between particles, ~ n _1 , both decrease as 
./V -1 so that their ratio is independent of particle number. 

The fact that the chaotic frozen- TV orbits appear to be- 
come "more nearly regular" as TV increases is consistent 
with recent claims that the "scale" associated with TV- 
body chaos decreases with increasing TV. Specifically, by 
comparing trajectories associated with two nearby initial 
conditions evolved in the same frozen- TV potential, Val- 
luri and Merritt |Q found that, when scaled in terms 
of i? sys , the size of the system, the typical separation 
R sa t on which the initial exponential divergence satu- 
rates decreases with increasing particle number, so that 
Rsat/Rsys is a decreasing function of TV. 

That the rate of divergence in the nonlinear regime 
slows more and more for larger TV can be quantified 
by tracking the actual evolution of two orbits generated 
from nearby initial conditions and determining the time 
required before their separation becomes 'macroscopic' 
The result of such an investigation is illustrated in FIG- 
URE 12, which was generated once again from ensembles 
of 20 frozen-TV orbits all evolved with e = 1CP 5 . In each 
case, the unperturbed orbits were identical to those used 
to generate FIGURE 1; the perturbed orbits involved 
changing the initial value of x by an amount Sx = 10~ 6 . 
FIGURE 12 exhibits as a function of TV the mean time r 
required before the separation 

6r={6x 2 +5y 2 + 5z 2 ) 1 / 2 (5.3) 

had achieved the value Sr = 1. (For this initial condi- 
tion, the average value of r associated with the smooth 
characteristic was R s » 1.83.) The error bars were de- 
rived by considering the first and second ten orbits in 
the ensemble separately. Because individual orbits di- 
verge at vastly different rates, the dispersion associated 
with a 20 orbit ensemble is much larger than reflected 
by these error bars. It is clear that r increases system- 
atically with increasing TV, although considerably more 
slowly than with the TV 1 / 2 dependence observed for the 
divergence time scales tc and tj. 

When viewed in terms of collisionless invariants like an- 
gular momentum, the divergence of frozen- TV orbits from 
smooth characteristics with the same initial condition is 
well approximated as a diffusion process, in which AJ 
grows as (t/tj) 1 / 2 and where, for fixed tr>, the divergence 
time scale tj varies at least approximately as TV 1 / 2 . This 
reinforces the conventional wisdom M| that discreteness 
effects may be modeled as white, or nearly white, Gaus- 
sian noise in the context of a Langevin or Fokker-Planck 
description. It might, therefore, seem somewhat surpris- 
ing that, although the divergence time scale to in config- 
uration or velocity space again scales as TV 1 / 2 , the quan- 
tities Ar and Av grow linearly in time, rather than as 
t 1 / 2 . 

In this regard, it is significant that if the smooth Plum- 
mer potential be replaced by the smooth potential asso- 
ciated with a constant density configuration, the linear 



growth exhibited Ar and Av is in fact replaced by the 
'expected' diffusive behaviour. In this case Ar and Av 
both grow as t 1 / 2 , and, when expressed in units of the dy- 
namical time to, the growth time to is somewhat longer, 
corresponding more nearly to the time scale tj associated 
with AJ. 

This suggests strongly that the behaviour of Ar and 
Av observed for the Plummer potential is associated with 
linear phase mixing. Because of finite number statistics, 
the same initial condition (r ,v ) in different frozen- TV 
realisations of a Plummer potential will correspond to 
somewhat different energies, the values of which are con- 
served in the subsequent evolution. However, even ne- 
glecting discreteness effects, initially proximate orbits in 
a generic integrable potential will, if their energies be un- 
equal, tend to diverge linearly. For example, two orbits 
evolved in a smooth Plummer potential with the same 
initial r but slightly different values of v and, hence, 
slightly different energies, will oscillate with somewhat 
different frequencies and, as a result, exhibit an overall 
linear divergence. If, however, the orbits are evolved in- 
stead in the potential associated with a constant density 
distribution, this is no longer true. A constant density 
sphere corresponds to a harmonic potential, where all 
orbits have the same unperturbed frequencies; and, for 
this reason, orbits in the smooth potential with slightly 
different energies will not exhibit such a systematic di- 
vergence. 

The fact that frozen- TV orbits look "more nearly regu- 
lar" for large TV suggests that the chaos associated with 
discreteness effects in the TV-body problem should be 
viewed very differently from the chaos associated with 
a bulk nonintegrable potential. When evolved into the 
future, two nearby chaotic initial conditions in such a 
potential tend to diverge exponentially until they are 
separated by a distance comparable to the size of the 
easily accessible (i.e., not significantly impeded by the 
Arnold web) connected phase space region to which the 
orbits are confined, a region which tends, typically, to 
be macroscopic. By contrast, the scale associated with 
chaos induced by discreteness effects in the TV-body prob- 
lem is distinctly microscopic, at least for comparatively 
large TV. It would appear that any single orbit with fixed 
energy can access a phase space region which is in fact 
very large; but the chaos which it experiences is a super- 
position of short range effects with characteristic scale 

Rsys- 
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FIG. 1. Mean short time Lyapunov exponent (x) as a 
function of softening parameter e for N = 10 5 (solid line), 
N = 10 4 ' 5 (dotted), N = 10 4 (dashed), N = W 3S 
(dot-dashed), N = 10 3 (triple-dot dashed), and N = 10 2 ' 5 
(broad dashed). The integrations were all performed for a 
single 'typical' initial condition. 
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FIG. 2. Mean short time Lyapunov exponent (x) as a func- 
tion of particle number N for e = 10~ 5 (solid line), e = 10~ 4 
(dotted), e = 10~ 3 (thin-dashed), e = 10~ 2 (dot-dashed), 
and e = 10 _1 (triple-dot dashed), all computed for the initial 
condition used to generate FIG. 1. The short time {\) for 
a different initial condition corresponding to a smooth radial 
orbit, again evolved with e = 10~ 5 , is indicated by the curve 
with thick dashes. 




FIG. 3. log 10 <7 x (Ai) as a function of log 10 At for three sets 
of simulations: N = 31623 and e — 0.0001 (solid curve), 
N = 316 and e = 0.0001 (dashed curve), and N = 316 and 
e = 0.1 (dot-dashed curve), all computed for the initial condi- 
tion used to generate FIG. 1. The thick solid line has a slope 
corresponding to p — 0.4. 



FIG. 5. (a) The radial coordinate r s (t) in the smooth po- 
tential (thin curve) and the mean trajectory (r(t)} (thick 
curve) derived from 20 frozen-JV simulations with N = 316 
and e = 10~ 4 , performed for an initial condition correspond- 
ing in the smooth potential to a purely radial orbit, (b) The 
same for N = 1000. (c) N = 3162. (d) N = 10000. (e) 
N = 31623. (f) N = 100000. 
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FIG. 4. (a) The trajectory x S m(t) in the smooth poten- 
tial (thin curve) and the mean trajectory (x(t)) (thick curve) 
derived from 20 frozen-iV simulations with iV = 316 and 
e = 10 -4 , performed for the initial condition used to gen- 
erate FIG 1. (b) The same for N = 1000. (c) N = 3162. (d) 
N = 10000. (e) N = 31623. (f) N = 100000. 
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FIG. 6. (a) The quantity Ar/V2R4 for frozen- N simula- 
tions with N = 1000 and e = 10 -4 . (b) The same data 
subjected to boxcar averaging over an interval t = tn- (c) 
and (d) The same for N = 100000. 
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FIG. 7. Best fit values of the time scale tc{N) associated 
with the divergence of smooth and frozen- N orbits for two dif- 
ferent initial conditions: Ar/R s — Av/V 3 = t/ta- The dashed 
line has slope 1/2, corresponding to an N 1 ' 2 dependence. 
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FIG. 8. Best fit values of the time scale tj(N) associated 
with changes in angular momentum for frozen-iV orbits for 
two different initial conditions: AJ 2 /J 2 = t/tj. The dashed 
line has slope 1/2, corresponding to an N 1 ^ 2 dependence. 



FIG. 9. The x-y projection of representative frozen-iV or- 
bits generated from the same initial condition, evolved for 
t = 25t D with e = 1CT 5 . (a) N = 316. (b) N = 1000. (c) 
N = 3163. (d) N = 10000. (e) N = 31623. (f) N = 100000. 
(g) N = 316228. (h) The x-y projection of the same initial 
condition evolved in the smooth potential 
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FIG. 12. The mean time required for two frozen- TV orbits 
separated initially by a distance Sr = 10~ 6 to achieve a macro- 
scopic separation 5r = 1. 



FIG. 10. (a) The Fourier transformed |a;(w)| for one 
frozen-TV integration of the initial condition used to gener- 
ate FIGURE 4, evolved with e = 1(T 5 and TV = 316. (b) 
The same for TV = 1000. (c) TV = 3162. (d) TV = 10000. (e) 
TV = 31623. (f) TV = 100000. (g) TV = 316228. (h) |a;(u;)| for 
a characteristic in the smooth potential with the same initial 
condition, with data recorded at the same intervals for the 
same total integration time. 



3.0 




FIG. 11. Two probes of the complexity of frozen- TV or- 
bits for an ensemble of orbits with the same initial condition 
evolved with e = 10~ 5 . The solid curve exhibits /o.i, the num- 
ber of frequencies which have power equal to at least 10% of 
the power in the peak frequencies. The dashed curve exhibits 
fco.95, the number of frequencies required to capture 95% of 
the total power. The vertical lines show /o.i and fco.95 for 
a smooth characteristic generated identically from the same 
initial condition, thus exhibiting the intrinsic limitations as- 
sociated with the discrete time series of data points. 



